home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / ode_qc_forward.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-17  |  4.8 KB  |  139 lines
  1. /*
  2. ### driver for quality controlled integrator ###
  3.  
  4. -----------------------------------------------------------------------------
  5. Note:    - Integration is done in standard coods, that is, the coordinate system defining
  6.     the original vector field.
  7.     - Plotting is done in window coordinates, the coord system as defined in
  8.     window_size panel 
  9. -----------------------------------------------------------------------------
  10. - int_driver:    0 no sense of absolute time, integration from i=0 to i_max
  11.         1 knows absolute time, integration from int_start to int_end
  12.         in int_nsteps;
  13.         2 knows absolute time, integration from int_start to int_end
  14.         with variable steps but error tolerance less than int_eps 
  15.  
  16. - draw_all:    0 compute section by sampling at uniformly spaced
  17.           discrete time intervals (int_period)
  18. -----------------------------------------------------------------------------
  19. - kount: counts the number of computed and displayed data (to be used for
  20. checking transients)
  21.   BS case is not fully implemented since the number of points computed
  22.   are usually small. The inclusion of intermediate points can be turned on
  23.   by setting int_option=1 but these points are not accurate)
  24. -----------------------------------------------------------------------------
  25. THIS DRIVER DEALS ONLY WITH THE CASE int_driver=2
  26. SECTION COMPUTATION NOT IN A DIFFERENT WAY
  27. -----------------------------------------------------------------------------
  28. */
  29.  
  30. #include <math.h>
  31. #define TINY 1.0e-30
  32.  
  33. int kount;
  34.  
  35. int
  36. ode_qc_forward()
  37. {
  38.     int i,nstep,int_t_end;
  39.     double int_t,int_next_step,int_actual_step,int_step;
  40.     double *vfull,*y,*dydx,*dvector(),fabs();
  41.     extern int nok,nbad,var_dim,func_dim,param_dim,full_dim;
  42.     extern int model,polar_coord,stop,draw_all;
  43.     extern int nok,nbad,int_driver,int_algorithm,int_max_nstep;
  44.     extern double int_start,int_end,int_period,int_eps,int_guessed_step,int_min_step,*int_yscal;
  45.     extern double *win_var_i,*win_var_f,*param;
  46.     extern char string[];
  47.     extern int (*f_p)();
  48.  
  49.  
  50.     stop = 0;
  51.     /* memory allocations */
  52.     y = dvector(0,var_dim-1);
  53.     dydx = dvector(0,var_dim-1);
  54.     vfull= dvector(0,full_dim-1);
  55.  
  56.     kount = 0;
  57.     nok = 0;
  58.     nbad = 0;
  59.     int_step = (int_end > int_start) ? fabs(int_guessed_step) : -fabs(int_guessed_step);
  60.     int_t = int_start;
  61.  
  62.     /* intialize from starting window variables to standard variables */
  63.     from_window_variables(y,win_var_i,polar_coord);
  64.     to_full_variables(vfull,y,polar_coord);
  65.     if(draw_all){
  66.         (void) draw_record_orbit(vfull,0,0);
  67.     }
  68.  
  69. restart:
  70.     if(!draw_all)
  71.         int_t_end = int_t+int_period;
  72.     else
  73.         int_t_end = int_end;
  74.         
  75.     for(nstep = 1; nstep <= int_max_nstep; nstep++) {
  76.         (int) f_p(dydx,0,y,param,int_t,var_dim);
  77.         /* set the error scaling function */
  78.         for(i=0;i<var_dim;i++)
  79.             int_yscal[i] = fabs(y[i]) + fabs(dydx[i] * int_step) + TINY;
  80.         /* do not draw and record intermediate orbits in the case of section */ 
  81.         if(draw_all){
  82.             to_full_variables(vfull,y,polar_coord);
  83.             (void) draw_record_orbit(vfull,kount++,1);
  84.         }
  85.  
  86.         /* set the final integration step of each
  87.         integration segment to end at integer multiples 
  88.         of int_period after t_start */
  89.         if((int_t + int_step - int_t_end) * (int_t + int_step - int_start) > 0.0) int_step = int_t_end -int_t;
  90.  
  91.         /* choose either runge-kutta or Bulirsch-Stoer quality controlled algorithm */
  92.         if(int_algorithm == 0)
  93.             (void) rkqc(y,dydx,var_dim,&int_t,int_step,int_eps,int_yscal,&int_actual_step,&int_next_step);
  94.         else if(int_algorithm == 1)
  95.             (void) bsstep(y,dydx,var_dim,&int_t,int_step,int_eps,int_yscal,&int_actual_step,&int_next_step);
  96.         
  97.         /* increment the number of bad or good steps */
  98.         if(int_actual_step == int_step) ++nok; else ++nbad;
  99.  
  100.         /* Record the final data in each integration segment.
  101.            While in the middle of integration segments
  102.            restart the case of section. Otherwise exit. */
  103.         if((int_t- int_t_end) * (int_t_end - int_start) >= 0.0) {
  104.             to_full_variables(vfull,y,polar_coord);
  105.             (void) draw_record_orbit(vfull,kount++,1);
  106.             /* in case of section, continue until t < t(end) */
  107.             if(!draw_all && int_t < int_end)
  108.                 goto restart;
  109.             else
  110.                 goto done;
  111.         }
  112.  
  113.         /* if integration step is smaller than int_min_step exit */
  114.         if( fabs(int_next_step) <= int_min_step) {
  115.             system_mess_proc(1,"ode_qc_forward: Step size too small.Stop!");
  116.             goto done;
  117.         }
  118.         if(stop){
  119.             goto done;
  120.         }    
  121.         int_step = int_next_step;
  122.     }
  123.     system_mess_proc(1,"ode_qc_forward: Too many steps. Stop!");
  124.  
  125. done:
  126.     /* reset the final value of window variables */
  127.     to_window_variables(win_var_f,y,polar_coord);
  128.     if(int_algorithm==0)
  129.         printf("Runge-Kutta QC Integ: nok=%d,nbad=%d\n",nok,nbad);
  130.     else if(int_algorithm == 1)
  131.         printf("Bulirsch-Stoer QC Integ: nok=%d,nbad=%d\n",nok,nbad);
  132.     (void) free_dvector(dydx,0,var_dim-1);
  133.     (void) free_dvector(y,0,var_dim-1);
  134.     (void) free_dvector(vfull,0,full_dim-1);
  135.     /* should be changed to return "time" instead ; return(int_t) */
  136.     return(kount);
  137. }
  138.  
  139.